library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(DT)
library(ggpubr)
library(cowplot)
source("scripts/utils/utils.R")
modernacolors()## <environment: R_GlobalEnv>
# load ELISpot data
temp <- read_excel("processed_data/ELISpot.xlsx")
dp <- temp[!is.na(temp$counts), ] %>%
filter(`low cell number (<0.15M)`== "no",
day != "143")%>%
mutate(counts = counts + 0.001)
d <-
dp %>%
filter(interval != "PBS")%>%
mutate(animal_id = factor(animal_id),
interval = factor(interval),
day = factor(day))dp %>%
filter(interval != "PBS") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
ggplot(aes(x = interval, y = counts)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "red",
width = .5,
geom = "errorbar",
size = 1
) +
facet_wrap(tissue ~ coating, scales = "free") +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "antigen specific \n ASC counts/million cells",
color = "",
linetype = "",
title = "Antigen specific spleen and BM results"
) +
theme(
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
legend.position = "bottom"
)#-------------------------------------------------------------------------------
# Fit the model for spleen data
#-------------------------------------------------------------------------------
# boxcox function to fit boxcox regression
boxcox_fit <- function(dat) {
require(MASS)
if (length(unique(dat$day)) == 1) {
out <- boxcox(counts ~ interval, data = dat)
optimal <- out$x[which.max(out$y)]
bctran <- make.tran("boxcox", optimal)
lm_fit <- with(bctran,
lm(linkfun(counts) ~ interval, data = dat))
} else{
out <- boxcox(counts ~ interval * day, data = dat)
optimal <- out$x[which.max(out$y)]
bctran <- make.tran("boxcox", optimal)
lm_fit <- with(bctran,
lm(linkfun(counts) ~ interval * day, data = dat))
}
return(list(lm_fit = lm_fit, lambda = optimal))
}
# box cox model fit for spleen data
d_tissue <- subset(d, tissue == "spleen") %>%
filter(!(coating == "NTD" & interval == "0"))
boxcox_fits <-
d_tissue %>%
split(.$coating) %>%
map( ~ boxcox_fit(.x))lm_fits <-
purrr::transpose(boxcox_fits)$lm_fit
# obtain optimal lambda for boxcox transformation
lambdas <-
purrr::transpose(boxcox_fits)$lambda %>% unlist()
assumption_plots <-
imap(
lm_fits ,
~ assumption_check(.x, .x$model$`linkfun(counts)`) %>%
plot_grid(
ggdraw() + draw_label(
paste0("Spleen ", .y, "-specific model residual diagnostics "),
x = 0.03,
hjust = 0
),
.,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
)
# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
emmeans(lm_grid, specs = ~ interval, type = "response") %>%
as_tibble()ggarrange(plotlist = assumption_plots)#-------------------------------------------------------------------------------
# Figure 4A. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
# fold change figure
foldchange1 <-
emmt_dat %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(facet = "1 week post-boost") %>%
ggplot(aes(x = interval , y = response, color = interval)) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.5),
size = 1
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
#ggtitle("S2P specific ASCs") +
facet_grid( ~ facet) +
labs(x = "", y = "Estimated counts/million cells") +
geom_hline(yintercept = emmt_dat$lower.CL[c(6, 7)],
color = "darkcyan",
lty = "dashed") +
theme(
axis.text = element_text(size = 12),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
# boxplot
box1 <-
dp %>%
filter(interval != "PBS", tissue == "spleen", coating == "S2P") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(facet = "1 week post-boost") %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
ggplot(aes(x = interval, y = counts, color = interval)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "grey",
width = .5,
geom = "errorbar",
size = 1
) +
facet_grid( ~ facet) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "Counts/million cells",
color = "",
linetype = "",
title = "S2P specific ASCs"
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme(
axis.text = element_text(size = 12),
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
legend.position = "None",
)#-------------------------------------------------------------------------------
# Supplementary Table 3
#-------------------------------------------------------------------------------
emm_elispot <-
emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
specs = ~ interval,
type = "response")
set.seed(100)
sum_elispot <-
contrast(
regrid(emm_elispot, transform = "log10"),
"pairwise",
infer = TRUE,
adjust = "mvt",
ratio = FALSE,
level = .95
) %>% as_tibble() %>%
bind_rows()
ELISpot_P_ST3 <-
sum_elispot %>%
as_tibble() %>%
dplyr::select(contrast, estimate, p.value) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value"
) %>%
mutate(`Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
`Fold change`= round(`Fold change`, digits = 6))%>%
filter(`Adjusted P-value`< 0.05)
ELISpot_P_ST3 %>%
DT::datatable(caption = "Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing interval") # Supplementary Table 3 visualization
t3v <-
sum_elispot %>%
as_tibble() %>%
separate(contrast, c("a", "b"), sep = " - ") %>%
mutate(a = str_remove_all(a , "interval"),
b = str_remove_all(b , "interval")) %>%
mutate(
a = case_when(
str_detect(a, "1") ~ paste(a, "week"),
str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
TRUE ~ paste(a, "weeks")
),
b = case_when(
str_detect(b, "1") ~ paste(b, "week"),
str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
TRUE ~ paste(b, "weeks")
)
) %>%
mutate(contrast = paste(a, "/", b)) %>%
mutate(day = "1 week post-boost)") %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
) %>%
mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
#position = position_dodge(width = 2),
size = 1) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
facet_grid( ~ day) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(title = "S2-P-specific ASCs",
y = "Estimated Fold-Change",
x = "Prime-boost dosing interval") +
coord_flip() +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 12),
strip.text = element_text(size = 10),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)#-------------------------------------------------------------------------------
# Fit the model for BM data
#-------------------------------------------------------------------------------
# BM
d_bm <- subset(d, tissue == "BM")
boxcox_fits <-
d_bm %>%
split(.$coating) %>%
map(~ boxcox_fit(.x))lm_fits <-
purrr::transpose(boxcox_fits)$lm_fit
# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
emmeans(lm_grid, specs = ~ interval | day, type = "response") %>%
as_tibble()assumption_plots <-
imap(
lm_fits ,
~ assumption_check(.x, .x$model$`linkfun(counts)`) %>%
plot_grid(
ggdraw() + draw_label(
paste0("BM ", .y, "-specific model residual diagnostics "),
x = 0.03,
hjust = 0
),
.,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
)
ggarrange(plotlist = assumption_plots)Each dot in the bottom panel of Figure 4 corresponds to the estimated mean for the prime-boost dosing interval. Error bars are 95% CIs on the estimated mean.
#-------------------------------------------------------------------------------
# Figure 4B. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
foldchange2 <-
emmt_dat %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(
factor(day),
"4 week post-boost" = "87",
"24 week post-boost" = "226"
)) %>%
mutate(day = factor(day, c("4 week post-boost" , "24 week post-boost"))) %>%
ggplot(aes(x = interval , y = response, color = interval)) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.5),
size = 1
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
labs(x = "", y = "Estimated counts/million cells") +
geom_hline(yintercept = emmt_dat$lower.CL[c(6)],
color = "darkcyan",
lty = "dashed") +
geom_hline(yintercept = emmt_dat$lower.CL[c(13)],
color = "darkcyan",
lty = "dashed") +
facet_grid(~ day) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme(
axis.text = element_text(size = 12),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_blank(),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
box2 <-
dp %>%
filter(interval != "PBS", tissue == "BM", coating == "S2P") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(
factor(day),
"4 week post-boost" = "87",
"24 week post-boost" = "226"
)) %>%
mutate(day = factor(day,
levels = c("4 week post-boost", "24 week post-boost"))) %>%
ggplot(aes(x = interval, y = counts, color = interval)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "grey",
width = .5,
geom = "errorbar",
size = 1
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "Counts/million cells",
color = "",
linetype = "",
title = "S2-P-specific LLPCs"
) +
facet_grid(~ day) +
theme(
axis.text = element_text(size = 12),
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold")
)
cowplot::plot_grid(
box1,
box2,
foldchange1,
foldchange2,
labels = c("A", "C", "B", "D"),
align = 'hv',
rel_widths = c(1, 2),
rel_heights = c(1, 1)
)#-------------------------------------------------------------------------------
# Supplementary Table 4
#-------------------------------------------------------------------------------
emm_elispot_BM <-
emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
specs = ~ interval |
day,
type = "response")
set.seed(100)
emm_elispot_BM <-
contrast(
regrid(emm_elispot_BM , transform = "log10"),
"pairwise",
infer = TRUE,
adjust = "mvt",
ratio = FALSE,
level = .95
) %>% as_tibble() %>%
bind_rows()
ELISpot_P_ST4 <-
emm_elispot_BM %>%
as_tibble() %>%
dplyr::select(day, contrast, estimate, p.value) %>%
mutate(day = as.numeric(as.character(day))) %>%
arrange(day) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
"Day" = "day"
) %>%
mutate(`Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
`Fold change`= round(`Fold change`, digits = 6))%>%
filter(`Adjusted P-value`< 0.05)
ELISpot_P_ST4 %>%
DT::datatable(caption = "Statistical comparisons of S2-P–specific LLPCs 4 weeks following
dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273
(10 μg) dosing intervals") Each dot corresponds to the estimated fold change for the contrasted prime-boost dosing interval. Error bars are multivariate t adjusted 95% CIs on the estimated contrast.
t4v <-
emm_elispot_BM %>%
as_tibble() %>%
separate(contrast, c("a", "b"), sep = " - ") %>%
mutate(a = str_remove_all(a , "interval"),
b = str_remove_all(b , "interval")) %>%
mutate(
a = case_when(
str_detect(a, "1") ~ paste(a, "week"),
str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
TRUE ~ paste(a, "weeks")
),
b = case_when(
str_detect(b, "1") ~ paste(b, "week"),
str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
TRUE ~ paste(b, "weeks")
)
) %>%
mutate(contrast = paste(a, "/", b)) %>%
mutate(day = factor(day, levels = c(87, 226))) %>%
mutate(day = forcats::fct_recode(
day,
"4 week post-boost" = "87",
"24 week post-boost" = "226",
)) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
) %>%
mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
#position = position_dodge(width = 2),
size = 1) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
facet_grid(~ day) +
labs(title = "S2-P-specific LLPC",
# subtitle = "4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227)
# between mRNA-1273 (10 μg) dosing intervals",
y = "Estimated Fold-Change",
x = "Prime-boost dosing interval") +
coord_flip() +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 12),
strip.text = element_text(size = 10),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
cowplot::plot_grid(
plotlist = list(t3v, t4v),
labels = c("A", "B"),
nrow = 1,
ncol = 2,
rel_widths = c(0.5, 1)
)sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MASS_7.3-58.1 cowplot_1.1.1 ggpubr_0.4.0 DT_0.24
## [5] lme4_1.1-30 Matrix_1.5-1 ggeffects_1.1.3 emmeans_1.8.0
## [9] mgcv_1.8-40 nlme_3.1-159 readxl_1.4.1 forcats_0.5.2
## [13] stringr_1.4.1 dplyr_1.0.10 purrr_0.3.4 readr_2.1.2
## [17] tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 lubridate_1.8.0 httr_1.4.4
## [4] tools_4.2.1 backports_1.4.1 bslib_0.4.0
## [7] utf8_1.2.2 R6_2.5.1 DBI_1.1.3
## [10] colorspace_2.0-3 withr_2.5.0 tidyselect_1.1.2
## [13] compiler_4.2.1 cli_3.3.0 rvest_1.0.3
## [16] xml2_1.3.3 labeling_0.4.2 sass_0.4.2
## [19] scales_1.2.1 mvtnorm_1.1-3 digest_0.6.29
## [22] minqa_1.2.4 rmarkdown_2.16 pkgconfig_2.0.3
## [25] htmltools_0.5.3 highr_0.9 dbplyr_2.2.1
## [28] fastmap_1.1.0 htmlwidgets_1.5.4 rlang_1.0.5
## [31] rstudioapi_0.14 jquerylib_0.1.4 generics_0.1.3
## [34] farver_2.1.1 jsonlite_1.8.0 crosstalk_1.2.0
## [37] car_3.1-0 googlesheets4_1.0.1 magrittr_2.0.3
## [40] Rcpp_1.0.9 munsell_0.5.0 fansi_1.0.3
## [43] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.8
## [46] yaml_2.3.5 carData_3.0-5 grid_4.2.1
## [49] crayon_1.5.1 lattice_0.20-45 haven_2.5.1
## [52] splines_4.2.1 hms_1.1.2 knitr_1.40
## [55] pillar_1.8.1 boot_1.3-28 estimability_1.4.1
## [58] ggsignif_0.6.3 reprex_2.0.2 glue_1.6.2
## [61] evaluate_0.16 modelr_0.1.9 vctrs_0.4.1
## [64] nloptr_2.0.3 tzdb_0.3.0 cellranger_1.1.0
## [67] gtable_0.3.1 assertthat_0.2.1 cachem_1.0.6
## [70] xfun_0.32 xtable_1.8-4 broom_1.0.1
## [73] coda_0.19-4 rstatix_0.7.0 googledrive_2.0.0
## [76] gargle_1.2.0 writexl_1.4.0 ellipsis_0.3.2